Estimated change in pneumonia in Latin American countries, subnational analyses

Daniel Weinberger

2019-08-08


Installing the package

Uncomment and run the following 3 lines of code to install the package. You will be prompted to manually select a number to indicate whether to update dependencies and which ones.

#install.packages('devtools') 
#library(devtools) 
#devtools::install_github('https://github.com/weinbergerlab/InterventionEvaluatR') 

Analysis goal

The goal for this analysis is to quantify changes in the rates of hospitalization due to pneumonia in different regions of Latin American countries


View dataset

Let’s load the data, which are included with the package

    ds<-read.csv('./prelog_Brazil_state_processed_data.csv')
  head(ds)
#>   age_group       date J12_18 P00_99 P05_07 a10_b99_nopneumo ach_noj
#> 1   02_1_11 2003-01-01     69     NA     NA               10     221
#> 2   02_1_11 2003-02-01     56     NA     NA               15     198
#> 3   02_1_11 2003-03-01     83     NA     NA               17     198
#> 4   02_1_11 2003-04-01     79     NA     NA               13     165
#> 5   02_1_11 2003-05-01    110     NA     NA               22     214
#> 6   02_1_11 2003-06-01    116     NA     NA               18     213
#>   K00_99 cj20_J22 A41 E00_99 N00_99 Q00_99 Z00_99 G00_99_SY E40_46 L00_99
#> 1     10       NA  NA     NA     NA     NA     NA        NA     NA     NA
#> 2     12       NA  NA     NA     NA     NA     NA        NA     NA     NA
#> 3     12       NA  NA     NA     NA     NA     NA        NA     NA     NA
#> 4     10       NA  NA     NA     NA     NA     NA        NA     NA     NA
#> 5     16       NA  NA     NA     NA     NA     NA        NA     NA     NA
#> 6     13       NA  NA     NA     NA     NA     NA        NA     NA     NA
#>   S00_T99 I00_99 N39 C00_D48 D50_89 H00_99_SY M00_99 E10_14 A39 B34 B20_24
#> 1      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#> 2      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#> 3      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#> 4      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#> 5      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#> 6      NA     NA  NA      NA     NA        NA     NA     NA  NA  NA     NA
#>   K35 I60_64 K80 V00_Y99 A18 B96
#> 1  NA     NA  NA      NA  NA  NA
#> 2  NA     NA  NA      NA  NA  NA
#> 3  NA     NA  NA      NA  NA  NA
#> 4  NA     NA  NA      NA  NA  NA
#> 5  NA     NA  NA      NA  NA  NA
#> 6  NA     NA  NA      NA  NA  NA
  #ds$G00_99_SY <-NULL #Exclude meningitis as a control 
ds<-ds[substr(ds$age_group,1,2) %in% c('09'),]

Ensure your date variable is in an R date format. If your variable is in a character or factor format, you need to tell R the format. – %m is a 2 digit month; %b is a month abbreviation (ie Jan, Feb) – %d is a 2 digit day (0-31), – %Y is a 4 digit year (e.g. 2011), %y is a 2 digit year (e.g. 11).
These codes are separated by a dash or slash or space. Modify the tryFormats script below if needed to match the format in your dataset

ds$date<-as.Date(ds$date, tryFormats=c('%Y-%m-%d',
                                                    '%m-%d-%Y',
                                                    '%m/%d/%Y',
                                                    '%Y/%m/%d',
                                                    '%d/%m/%Y'
                                                    ) )

Sort data by age group and month

    ds<-ds[order(ds$age_group, ds$date),] #Sort data by age group and month

Set parameters for analysis

Here we need to set a few parameters. We Use the evaluatr.init() function to specify the name of the dataset, the date at which the vaccine is introduced, the date at which we want to begin evaluating the vaccine (typically 1-2 year after vaccine introduction). We also provide some information on the dataset, sch as whether the data are monthly or quarterly (n_seasons), the variable names for the grouping variable, the date variable, the outcome variable, and the denominator variable (if any).


analysis <- evaluatr.init(
  country = "Brazil_state", data = ds,
  post_period_start = "2010-03-01", #First 'post-intervention' month is Jan 2012
  eval_period_start = "2011-03-01", #We ignore first 12 month of data to allow for vaccine ramp up
  eval_period_end = "2014-12-01", #The evaluation period lasts 2 years
  n_seasons = 12, #This is monthly data, so select 12
  year_def = "cal_year", # we are in southern hemisphere, so aggregate results by calendar year (Jan-Dec)
  group_name = "age_group",  #Strata categry name
  date_name = "date", #Date variable name
  outcome_name = "J12_18", #Outcome variable name
  denom_name = "ach_noj", #Denominator variable name
   set.burnN=10000,
  set.sampleN=5000
)
set.seed(1)

Run the main analysis

Save the results in object ‘impact_results’

Check convergence

n.iter<-ncol(impact_results$full$rr_iter)
thin=1
trace1<-impact_results$full$rr_iter[,seq(from=1, to=n.iter, by=thin)] #thin
for(i in 1: nrow(trace1)){
  geweke.p<- pnorm(abs(geweke.diag(mcmc(trace1[i,]))$z),lower.tail=FALSE)*2
  print(geweke.p)
if(geweke.p>0.05){con.stat<-'Model converged'
  }else{
    con.stat<-'Not converged'
  } 
   plot(trace1[i,], type='l', main=paste0(con.stat,' ',colnames(trace1)[i]), bty='l', ylim=c(0.2,2))
}
#>      var1 
#> 0.7213792

#>      var1 
#> 0.7820469

#>      var1 
#> 0.2596145

#>       var1 
#> 0.07196071

#>      var1 
#> 0.5383267

#>      var1 
#> 0.5402461

#>       var1 
#> 0.08662169

#>      var1 
#> 0.7301025

#>       var1 
#> 0.02421471

#>      var1 
#> 0.3326364

#>       var1 
#> 0.03235626

#>        var1 
#> 0.007394855

#>      var1 
#> 0.1087894

#>       var1 
#> 0.07885182

#>      var1 
#> 0.7732755

#>     var1 
#> 0.383947

#>        var1 
#> 9.91779e-05

#>      var1 
#> 0.1807919

#>      var1 
#> 0.2018847

#>       var1 
#> 0.06964387

#>       var1 
#> 0.09856711

#>        var1 
#> 0.001022382

#>      var1 
#> 0.5661093

#>      var1 
#> 0.1605447

#>      var1 
#> 0.9276389

#>      var1 
#> 0.5848874

#>      var1 
#> 0.5304436

#>      var1 
#> 0.6631171

#>      var1 
#> 0.5856675

#>       var1 
#> 3.6141e-08

#>       var1 
#> 0.02549698

#>        var1 
#> 0.002741437

#>      var1 
#> 0.3239681

Results

Compare estimates from different models

This shows the estimated rate ratio and 95% credible intervals from a synthetic controls analysis; a time-trend analysis where we used the specified denominator (all non-respiratory deaths) to adjust the number of pneumonia deaths in each month and a linear trend for time; a classic interrupted time series analysis (segmented regression); and the STL+PCA approach, which smooths and combines the control variables prior to including them in the model.

results.table<- cbind.data.frame(
  #impact_results$best$rr_mean_intervals, 
  impact_results$full$rr_mean_intervals, 
  impact_results$time$rr_mean_intervals, 
  #impact_results$time_no_offset$rr_mean_intervals, 
  impact_results$its$rr_mean_intervals, 
  impact_results$pca$rr_mean_intervals)

  table<-xtable(results.table)
  htmlTable(table)
Synthetic controls Estimate (95% CI) Time trend Estimate (95% CI) Classic ITS (95% CI) STL+PCA Estimate (95% CI)
09_1_11 0.88 (0.72, 1.05) 0.91 (0.68, 1.21) 1.02 (0.76, 1.37) 0.6 (0.43, 0.84)
09_1_12 0.71 (0.58, 0.96) 0.96 (0.68, 1.32) 0.96 (0.72, 1.29) 1.57 (0.97, 2.37)
09_1_13 1.08 (0.95, 1.2) 1.19 (0.88, 1.56) 0.94 (0.72, 1.24) 0.86 (0.74, 1)
09_1_14 1.24 (0.94, 1.6) 0.8 (0.56, 1.11) 0.57 (0.42, 0.76) 0.67 (0.5, 0.88)
09_1_15 0.71 (0.65, 0.78) 0.98 (0.82, 1.15) 0.78 (0.67, 0.9) 0.77 (0.69, 0.86)
09_1_16 1.15 (0.91, 1.41) 0.8 (0.58, 1.1) 1.03 (0.76, 1.4) 1.06 (0.84, 1.34)
09_1_17 0.8 (0.56, 0.99) 0.54 (0.41, 0.69) 0.87 (0.68, 1.12) 0.76 (0.51, 1.14)
09_1_AA 0.87 (0.81, 0.94) 0.96 (0.83, 1.1) 0.87 (0.77, 0.99) 0.9 (0.81, 0.99)
09_2_21 0.9 (0.77, 1.3) 0.9 (0.65, 1.14) 1.03 (0.85, 1.25) 1 (0.84, 1.19)
09_2_22 0.69 (0.45, 0.93) 0.86 (0.65, 1.1) 0.84 (0.7, 1.01) 0.71 (0.62, 0.82)
09_2_23 0.73 (0.62, 0.84) 0.79 (0.6, 0.99) 0.77 (0.63, 0.93) 0.73 (0.63, 0.86)
09_2_24 0.54 (0.47, 0.61) 0.52 (0.43, 0.64) 0.64 (0.51, 0.79) 0.61 (0.51, 0.73)
09_2_25 0.78 (0.66, 0.87) 0.66 (0.52, 0.84) 0.78 (0.67, 0.91) 0.74 (0.62, 0.87)
09_2_26 0.86 (0.7, 0.99) 0.65 (0.47, 0.84) 0.78 (0.63, 0.96) 0.84 (0.71, 0.99)
09_2_27 0.8 (0.6, 1.02) 0.56 (0.43, 0.67) 0.72 (0.58, 0.91) 0.76 (0.64, 0.92)
09_2_28 0.47 (0.4, 0.6) 0.47 (0.36, 0.63) 0.9 (0.65, 1.25) 0.43 (0.37, 0.5)
09_2_29 0.57 (0.52, 0.74) 0.7 (0.63, 0.77) 0.81 (0.72, 0.91) 0.76 (0.71, 0.82)
09_2_AA 0.7 (0.65, 0.75) 0.75 (0.65, 0.84) 0.9 (0.8, 1.02) 0.84 (0.76, 0.92)
09_3_31 0.6 (0.53, 0.68) 0.71 (0.59, 0.83) 0.8 (0.69, 0.93) 0.74 (0.66, 0.83)
09_3_32 0.48 (0.42, 0.55) 0.59 (0.46, 0.72) 0.74 (0.59, 0.92) 0.76 (0.57, 0.99)
09_3_33 0.75 (0.64, 0.86) 0.59 (0.49, 0.71) 0.59 (0.5, 0.69) 0.66 (0.57, 0.77)
09_3_35 0.91 (0.85, 0.98) 1.11 (0.95, 1.3) 0.98 (0.86, 1.12) 0.87 (0.8, 0.95)
09_3_AA 0.81 (0.78, 0.85) 0.86 (0.72, 0.98) 0.85 (0.76, 0.95) 0.84 (0.74, 0.93)
09_4_41 0.86 (0.61, 1.02) 0.97 (0.83, 1.14) 0.99 (0.88, 1.12) 0.96 (0.86, 1.07)
09_4_42 0.88 (0.78, 0.98) 0.99 (0.85, 1.17) 1.04 (0.91, 1.2) 1.39 (1.11, 1.73)
09_4_43 0.85 (0.79, 0.92) 1.02 (0.9, 1.15) 0.95 (0.83, 1.08) 0.91 (0.83, 1)
09_4_AA 0.77 (0.69, 0.9) 1.01 (0.9, 1.14) 0.99 (0.9, 1.09) 0.94 (0.86, 1.02)
09_5_50 0.83 (0.73, 0.93) 0.76 (0.62, 0.93) 0.75 (0.62, 0.92) 0.84 (0.72, 0.98)
09_5_51 0.87 (0.75, 0.98) 0.97 (0.79, 1.16) 1.03 (0.87, 1.22) 1.13 (0.9, 1.38)
09_5_52 0.66 (0.58, 0.75) 0.8 (0.63, 1.05) 0.94 (0.76, 1.15) 0.91 (0.75, 1.11)
09_5_53 0.49 (0.42, 0.57) 0.67 (0.5, 0.86) 0.66 (0.54, 0.82) 0.53 (0.45, 0.61)
09_5_AA 0.7 (0.63, 0.78) 0.85 (0.71, 1) 0.88 (0.76, 1.03) 0.9 (0.75, 1.05)
09_A_AA 0.77 (0.74, 0.8) 0.85 (0.77, 0.93) 0.88 (0.81, 0.96) 0.85 (0.78, 0.91)

Cases averted

How many cases were prevented from the time of vaccine introduction to the last time point in each stratum (+/- 95% CrI)? You can modify the number ‘last.point’ to pull out the cumulative number of cases at any point in the time series

last.point<-dim(impact_results$best$cumsum_prevented)[1]
cum.prevented<-impact_results$best$cumsum_prevented[last.point,,]

Number of variables selected in SC analysis

The Synthetic controls analysis uses a Bayesian variable selection algorithm to weight the candidate covariates. In each MCMC iteration, it tests a different combination of variables. The model size indicates how many variables are selected in any given model. If <1 variable is selected on average, this indicates that no suitable control variables were identified. In this example 1-2 variables were selected on average in the 2-23m and 2-59 m age categories, while no controls were identified in the 24-59m age group (the average model size is <1 (0.44)).

model_size = data.frame(t(analysis$model_size))
htmlTable(model_size, align='c')
X09_1_11 X09_1_12 X09_1_13 X09_1_14 X09_1_15 X09_1_16 X09_1_17 X09_1_AA X09_2_21 X09_2_22 X09_2_23 X09_2_24 X09_2_25 X09_2_26 X09_2_27 X09_2_28 X09_2_29 X09_2_AA X09_3_31 X09_3_32 X09_3_33 X09_3_35 X09_3_AA X09_4_41 X09_4_42 X09_4_43 X09_4_AA X09_5_50 X09_5_51 X09_5_52 X09_5_53 X09_5_AA X09_A_AA
1 2.21 2.83 3.2 2.23 2.43 1.7 0.49 2.21 3.53 3.1 2.8 1.35 1.58 5.15 1.4 2.64 2.78 4.36 4.01 3.11 2.96 5.56 3.32 3.42 3.38 3.24 4.23 2.71 3.76 2.05 2.63 3.87 2.98

Inclusion Probabilities

Here we can see which variables received the most weight in the models for each strata. Values range from 0-1, with values closer to 1 indicating that the variable was included in a larger proportion of the variable combinations that were tested A value of 1 would indicate that the variable was always included in the model, a value of 0.48 would indicate the variable was included in 48% of the models that were tested.

htmlTable(incl_probs)
Group Greatest Inclusion Variable Greatest Inclusion Probability Second Greatest Inclusion Variable Second Greatest Inclusion Probability Third Greatest Inclusion Variable Third Greatest Inclusion Probability
1 09_1_11 cj20_J22 1 ach_noj 0.58 P00_99 0.26
2 09_1_12 cj20_J22 0.97 ach_noj 0.93 P00_99 0.4
3 09_1_13 K00_99 1 cj20_J22 1 N00_99 0.7
4 09_1_14 ach_noj 1 P00_99 0.93 P05_07 0.31
5 09_1_15 cj20_J22 1 K00_99 0.99 E40_46 0.19
6 09_1_16 ach_noj 0.6 a10_b99_nopneumo 0.59 P05_07 0.29
7 09_1_17 ach_noj 0.15 P00_99 0.13 P05_07 0.06
8 09_1_AA cj20_J22 1 N39 1 C00_D48 0.03
9 09_2_21 cj20_J22 1 ach_noj 0.7 K00_99 0.68
10 09_2_22 a10_b99_nopneumo 0.8 P05_07 0.78 P00_99 0.78
11 09_2_23 cj20_J22 0.73 P00_99 0.63 Q00_99 0.35
12 09_2_24 ach_noj 0.7 L00_99 0.23 a10_b99_nopneumo 0.1
13 09_2_25 cj20_J22 0.99 A41 0.31 a10_b99_nopneumo 0.09
14 09_2_26 cj20_J22 1 A41 1 a10_b99_nopneumo 0.73
15 09_2_27 E40_46 0.57 A41 0.18 E00_99 0.17
16 09_2_28 cj20_J22 1 ach_noj 0.94 P00_99 0.22
17 09_2_29 P00_99 0.81 a10_b99_nopneumo 0.67 E40_46 0.4
18 09_2_AA cj20_J22 1 A41 1 E40_46 1
19 09_3_31 cj20_J22 1 E00_99 0.82 N39 0.74
20 09_3_32 cj20_J22 1 A41 0.99 N00_99 0.64
21 09_3_33 cj20_J22 1 E40_46 0.67 E00_99 0.28
22 09_3_35 cj20_J22 1 D50_89 1 E00_99 0.77
23 09_3_AA cj20_J22 1 L00_99 1 E00_99 0.82
24 09_4_41 cj20_J22 1 P00_99 0.83 K00_99 0.72
25 09_4_42 ach_noj 1 P00_99 1 cj20_J22 1
26 09_4_43 cj20_J22 1 E00_99 1 N39 0.99
27 09_4_AA cj20_J22 1 P00_99 1 E00_99 0.94
28 09_5_50 cj20_J22 1 ach_noj 0.61 P00_99 0.48
29 09_5_51 cj20_J22 1 K00_99 1 P00_99 0.99
30 09_5_52 cj20_J22 1 K00_99 0.41 Q00_99 0.18
31 09_5_53 cj20_J22 1 A41 0.62 a10_b99_nopneumo 0.41
32 09_5_AA cj20_J22 1 K00_99 0.74 N39 0.73
33 09_A_AA cj20_J22 1 E40_46 1 M00_99 0.92

Generate and save the plots

Plot the results for 1 age group

First look at the results from the synthetic controls model for 1 age group.

This first plot shows a monthly time series, with observed, fitted, and counterfacual values. The observed number of deaths is shown in the black line. The fitted values for the pre-vaccine period are shown in the red dotted line, and the counterfactual estimate with its 95% credible interval is shown as the white dotted line and gray shaded area. if the black line is below the gray shaded area, this would indicate that obsrved cases are lower than expected based on changes in the control diseases in the post-vaccine period. If the controls appropriately adjust for underlying trends, then this would reflect an effect of the vaccine.

      plots$groups[[1]]$pred_full 

It is sometimes easier to look at the results if we aggregate the observed and expected values up to the annual time scale. Here the observed values are shown as black dots. When the black dots go below the gray shaded area, this indicates that the observed cases are lower than expected based on changes in the control diseases in the post-vaccine period. If the controls appropriately adjust for underlying trends, then this would reflect an effect of the vaccine. The vertical line indicates the date of vaccine introduction. The year when the vaccine is introduced is included to the right of the line, even if the vaccine was introduced part-way through the year. For instance, regardless of whether the the vaccine was introduced in January 2009 or November 2009, the observed dot for 2009 would be to the right of the line.

      plots$groups[[1]]$pred_full_agg 

Finally, we can look at the cumulative cases prevented. In this example, there have been 445 cases prevented (95%CrI: 58, 931) from the time of vaccine introduction to the last day month of the study period. This is calculated by takin the difference between the observed and fitted number of cases in each month, and summing them. If atleast 1 control disease is identified from the synthetic controls model, then the result here is drawn from that model, otherwise, it is drawn from the STL+PCA model.

      plots$groups[[1]]$cumsum_prevented 

Printing plots for all models and age groups

We instead might want to just print everything for all age groups and models. We can use the following code to do that

Plot Observed vs expected yearly time series

for (group in names(plots$groups)) {
      par(mfrow=c(4,1))
      print(plots$groups[[group]]$pred_full_agg )
      print(plots$groups[[group]]$pred_best_agg )
      print(plots$groups[[group]]$pred_time_agg )
      print(plots$groups[[group]]$pred_pca_agg )
}

Use HDI intervals instead

The package calcuates both an equal-tail 95% credible interval and a 95% highest posterior density interval (HPD or HDI interval). The latter might be more appropriate when the posterior distribution is skewed. The objects with HDI intervals are saved with an appendix of ‘HDI’ so –impact_results\(best\)pred_quantiles gives the estimate for equal tail intervals while –impact_results\(best\)pred_quantiles_HDI gives the estimate for the HDI interval –impact_results\(best\)cumsum_prevented_hdi gives the estimate of cumulative cases prevented using HPI –impact_results\(best\)ann_pred_HDI gives the annual predicted cases with HDI intervals –impact_results\(best\)rr_mean_HDI gives the overall rate ratio caclulated with HDI intervals – impact_results\(best\)log_rr_hdi gives the pointswise log-RR with HDI intervals Here we compare the version calculated with equal-tailed intervals with HPD intervals. In this example they give very similar coverage

print("Equal tailed intervals")
#> [1] "Equal tailed intervals"
round(impact_results$best$rr_mean[,c(2,1,3)],2)
#>         Point Estimate Lower CI Upper CI
#> 09_1_11           0.88     0.72     1.05
#> 09_1_12           0.71     0.58     0.96
#> 09_1_13           1.08     0.95     1.20
#> 09_1_14           1.24     0.94     1.60
#> 09_1_15           0.71     0.65     0.78
#> 09_1_16           1.15     0.91     1.41
#> 09_1_17           0.76     0.51     1.14
#> 09_1_AA           0.87     0.81     0.94
#> 09_2_21           0.90     0.77     1.30
#> 09_2_22           0.69     0.45     0.93
#> 09_2_23           0.73     0.62     0.84
#> 09_2_24           0.54     0.47     0.61
#> 09_2_25           0.78     0.66     0.87
#> 09_2_26           0.86     0.70     0.99
#> 09_2_27           0.80     0.60     1.02
#> 09_2_28           0.47     0.40     0.60
#> 09_2_29           0.57     0.52     0.74
#> 09_2_AA           0.70     0.65     0.75
#> 09_3_31           0.60     0.53     0.68
#> 09_3_32           0.48     0.42     0.55
#> 09_3_33           0.75     0.64     0.86
#> 09_3_35           0.91     0.85     0.98
#> 09_3_AA           0.81     0.78     0.85
#> 09_4_41           0.86     0.61     1.02
#> 09_4_42           0.88     0.78     0.98
#> 09_4_43           0.85     0.79     0.92
#> 09_4_AA           0.77     0.69     0.90
#> 09_5_50           0.83     0.73     0.93
#> 09_5_51           0.87     0.75     0.98
#> 09_5_52           0.66     0.58     0.75
#> 09_5_53           0.49     0.42     0.57
#> 09_5_AA           0.70     0.63     0.78
#> 09_A_AA           0.77     0.74     0.80

print("HPD intervals")
#> [1] "HPD intervals"
round(impact_results$best$rr_mean_hdi,2)
#>         Point Estimate lower upper
#> 09_1_11           0.88  0.71  1.04
#> 09_1_12           0.71  0.55  0.91
#> 09_1_13           1.08  0.96  1.21
#> 09_1_14           1.24  0.93  1.57
#> 09_1_15           0.71  0.65  0.78
#> 09_1_16           1.15  0.91  1.41
#> 09_1_17           0.76  0.49  1.10
#> 09_1_AA           0.87  0.81  0.93
#> 09_2_21           0.90  0.76  1.28
#> 09_2_22           0.69  0.43  0.89
#> 09_2_23           0.73  0.62  0.83
#> 09_2_24           0.54  0.47  0.61
#> 09_2_25           0.78  0.66  0.87
#> 09_2_26           0.86  0.70  0.99
#> 09_2_27           0.80  0.61  1.03
#> 09_2_28           0.47  0.39  0.59
#> 09_2_29           0.57  0.52  0.74
#> 09_2_AA           0.70  0.65  0.75
#> 09_3_31           0.60  0.52  0.67
#> 09_3_32           0.48  0.41  0.55
#> 09_3_33           0.75  0.63  0.84
#> 09_3_35           0.91  0.86  0.98
#> 09_3_AA           0.81  0.77  0.85
#> 09_4_41           0.86  0.61  1.01
#> 09_4_42           0.88  0.78  0.97
#> 09_4_43           0.85  0.80  0.92
#> 09_4_AA           0.77  0.69  0.89
#> 09_5_50           0.83  0.73  0.93
#> 09_5_51           0.87  0.76  0.99
#> 09_5_52           0.66  0.58  0.75
#> 09_5_53           0.49  0.42  0.57
#> 09_5_AA           0.70  0.63  0.78
#> 09_A_AA           0.77  0.74  0.80

Compare pointwise coverage of the 95% CrI calculated as HDI or equal-tailed. In this example they are nearly identical. The red line denote the HPD intervals, the gray line denotes the equal-tailed

matplot(impact_results$best$pred_quantiles_HDI[,c(2,3),1], type='l', col='gray', lty=c(2,1,2), bty='l', ylim=c(0, max(impact_results$best$pred_quantiles_HDI[,,1])), ylab='Count')

matplot(impact_results$best$pred_quantiles[,c(1,3),1], type='l', col='red', lty=c(2,1,2), bty='l', ylim=c(0, max(impact_results$best$pred_quantiles_HDI[,,1])), ylab='Count', add=T)
points(analysis$input_data[, analysis$outcome_name])